clear; close all; clc;
%% 绘制 Fig.4 特征值变化曲线
%  - (a) 中参数设置：b0=1.3, b1=0.1, b2 变化，d2=0 ；即 b0、b1、d2 为标量，b2 为向量
%  - (b) 中参数设置：b2=1, b1=0.1, b0 变化，d2=0 ；即 b2、b1、d2 为标量，b0 为向量
%  - (c) 与 (d) 分别为 b1 变化时，直接传入 cos(z')=+1 或 -1 的情形


figure;

%% (a) b0=1.3, b1=0.1, b2∈[-3, 3], d2=0  => cos(z')=cos(0)=1
subplot(2,2,1);
b0 = 1.3;
b1 = 0.1;
d2 = 0;
% 这里让 b2 为变化向量
b2_vals = linspace(-3, 3, 500);

[Re2_a, Im2_a, Re3_a, Im3_a] = eigs_1(b0, b1, b2_vals, d2);
yline(1, 'k--', 'LineWidth', 1.2);
yline(-1, 'k--', 'LineWidth', 1.2);
hold on;
plot(b2_vals, Re2_a, 'r','LineWidth',1.5);
plot(b2_vals, Re3_a, 'b','LineWidth',1.5);
plot(b2_vals, Im2_a, 'r--','LineWidth',1.2);
plot(b2_vals, Im3_a, 'b--','LineWidth',1.2);
hold off; grid on;
xlabel('b_2');
ylabel('\lambda');
legend('\lambda_2 (real)','\lambda_3 (real)',...
       '\lambda_2 (imag)','\lambda_3 (imag)','Location','best');
title('(a) b_0=1.3, b_1=0.1, b_2\in[-3,3], d_2=0');

%% (b) b2=1, b1=0.1, b0∈[0, 6], d2=0  => cos(z')=cos(0)=1
subplot(2,2,2);
b2 = 1;
b1 = 0.1;
d2 = 0;
% 让 b0 为变化向量
b0_vals = linspace(0, 6, 500);

[Re2_b, Im2_b, Re3_b, Im3_b] = eigs_1(b0_vals, b1, b2, d2);
yline(1, 'k--', 'LineWidth', 1.2);
yline(-1, 'k--', 'LineWidth', 1.2);
hold on;
plot(b0_vals, Re2_b, 'r','LineWidth',1.5);
plot(b0_vals, Re3_b, 'b','LineWidth',1.5);
plot(b0_vals, Im2_b, 'r--','LineWidth',1.2);
plot(b0_vals, Im3_b, 'b--','LineWidth',1.2);
hold off; grid on;
xlabel('b_0');
ylabel('\lambda');
legend('\lambda_2 (real)','\lambda_3 (real)',...
       '\lambda_2 (imag)','\lambda_3 (imag)','Location','best');
title('(b) b_0\in[0,6], b_1=0.1, b_2=1, d_2=0');

%% (c) b0=1.3, b2=1, b1∈[-1, 0.5], d2=2n(1-b1)*pi, 取 cos(z')=1
subplot(2,2,3);
b0 = 1.3;
b2 = 1;
b1_vals = linspace(-1, 0.5, 500);
% 此情况直接令 cos(z')=1
cosZprime = 1;

[Re2_c, Im2_c, Re3_c, Im3_c] = eigs_2(b0, b1_vals, b2, cosZprime);
yline(1, 'k--', 'LineWidth', 1.2);
yline(-1, 'k--', 'LineWidth', 1.2);
hold on;
plot(b1_vals, Re2_c, 'r','LineWidth',1.5);
plot(b1_vals, Re3_c, 'b','LineWidth',1.5);
plot(b1_vals, Im2_c, 'r--','LineWidth',1.2);
plot(b1_vals, Im3_c, 'b--','LineWidth',1.2);
hold off; grid on;
xlabel('b_1');
ylabel('\lambda');
legend('\lambda_2 (real)','\lambda_3 (real)',...
       '\lambda_2 (imag)','\lambda_3 (imag)','Location','best');
title('(c) b_0=1.3, b_1\in[-1,0.5], b_2=1, d_2=2n(1-b_1)\pi');

%% (d) b0=1.3, b2=1, b1∈[-1, 0.5], d2=(2n-1)(1-b1)*pi, 取 cos(z')=-1
subplot(2,2,4);
b0 = 1.3;
b2 = 1;
b1_vals = linspace(-1, 0.5, 500);
% 此情况直接令 cos(z')=-1
cosZprime = -1;

[Re2_d, Im2_d, Re3_d, Im3_d] = eigs_2(b0, b1_vals, b2, cosZprime);
yline(1, 'k--', 'LineWidth', 1.2);
yline(-1, 'k--', 'LineWidth', 1.2);
hold on;
plot(b1_vals, Re2_d, 'r','LineWidth',1.5);
plot(b1_vals, Re3_d, 'b','LineWidth',1.5);
plot(b1_vals, Im2_d, 'r--','LineWidth',1.2);
plot(b1_vals, Im3_d, 'b--','LineWidth',1.2);
hold off; grid on;
xlabel('b_1');
ylabel('\lambda');
legend('\lambda_2 (real)','\lambda_3 (real)',...
       '\lambda_2 (imag)','\lambda_3 (imag)','Location','best');
title('(d) b_0=1.3, b_1\in[-1,0.5], b_2=1, d_2=(2n-1)(1-b_1)\pi');

%% 辅助函数：eigs_1
% 若 b0 为向量，则 b2 必须为标量；若 b2 为向量，则 b0 必须为标量。
function [Re2_vec, Im2_vec, Re3_vec, Im3_vec] = eigs_1(param1, b1, param2, d2)
% 计算对应的 特征值 λ₂ 和 λ₃ 的实部和虚部
    % 判断哪个参数为向量
    if isscalar(param1) && ~isscalar(param2)
        % 此时 param1 为标量（对应 b0），param2 为向量（对应 b2）
        scalar_val = param1;
        vec = param2;
        varyFlag = 'b2';  % 表示 b2 变化
    elseif ~isscalar(param1) && isscalar(param2)
        % 此时 param1 为向量（对应 b0），param2 为标量（对应 b2）
        vec = param1;
        scalar_val = param2;
        varyFlag = 'b0';  % 表示 b0 变化
    else
        error('必须有且只有一个参数为向量，其它为标量。');
    end

    % zprime = d2/(1 - b1)（这里 b1, d2 均为标量）
    zprime = d2/(1 - b1);
    % 预分配结果
    n = length(vec);
    Re2_vec = zeros(1, n);
    Im2_vec = zeros(1, n);
    Re3_vec = zeros(1, n);
    Im3_vec = zeros(1, n);

    for idx = 1:n
        switch varyFlag
            case 'b2'
                % b0 为标量，b2 变化
                b0_val = scalar_val;
                b2_val = vec(idx);
            case 'b0'
                % b2 为标量，b0 变化
                b0_val = vec(idx);
                b2_val = scalar_val;
        end
        delta = 4 * b0_val * b2_val * cos(zprime) + b1^2;
        if delta >= 0  
            sqrtVal = sqrt(delta);
            lambda2 = (b1 + sqrtVal)/2;
            lambda3 = (b1 - sqrtVal)/2;
            Re2_vec(idx) = lambda2;
            Im2_vec(idx) = 0;
            Re3_vec(idx) = lambda3;
            Im3_vec(idx) = 0;
        else
            sqrtValImag = sqrt(-delta);
            Re2_vec(idx) = b1/2;
            Im2_vec(idx) = sqrtValImag/2;
            Re3_vec(idx) = b1/2;
            Im3_vec(idx) = -sqrtValImag/2;
        end
    end
end

%% 辅助函数：eigs_2
% 该函数用于 (c) 与 (d) 情况，直接传入 cos(z') 的值（+1 或 -1）。
% 此时 b0, b2, cosZp 均为标量，b1 为向量。
function [Re2_vec, Im2_vec, Re3_vec, Im3_vec] = eigs_2(b0, b1Range, b2, cosZp)
    n = length(b1Range);
    Re2_vec = zeros(1, n);
    Im2_vec = zeros(1, n);
    Re3_vec = zeros(1, n);
    Im3_vec = zeros(1, n);

    for idx = 1:n
        b1_val = b1Range(idx);
        val = 4 * b0 * b2 * cosZp + b1_val^2;
        if val >= 0
            sqrtVal = sqrt(val);
            lambda2 = (b1_val + sqrtVal)/2;
            lambda3 = (b1_val - sqrtVal)/2;
            Re2_vec(idx) = lambda2;
            Im2_vec(idx) = 0;
            Re3_vec(idx) = lambda3;
            Im3_vec(idx) = 0;
        else
            sqrtValImag = sqrt(-val);
            Re2_vec(idx) = b1_val/2;
            Im2_vec(idx) = sqrtValImag/2;
            Re3_vec(idx) = b1_val/2;
            Im3_vec(idx) = -sqrtValImag/2;
        end
    end
end
